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LOCAL MINIMIZATION ALGORITHMS FOR DYNAMIC 
PROGRAMMING EQUATIONS 

DANTE KALISE*, AXEL KRONERt, AND KARL KUNISCRt 

Abstract. The numerical realization of the dynamic programming principle for continuous-time 
optimal control leads to nonlinear Hamilton-Jacobi-Bellman equations which require the minimiza¬ 
tion of a nonlinear mapping over the set of admissible controls. This minimization is often performed 
by comparison over a finite number of elements of the control set. In this paper we demonstrate 
the importance of an accurate realization of these minimization problems and propose algorithms by 
which this can be achieved effectively. The considered class of equations includes nonsmooth control 
problems with -penalization which lead to sparse controls. 

Key words. dynamic programming, Hamilton-Jacobi-Bellman equations, semi-Lagrangian 
schemes, first order primal-dual methods, semismooth Newton methods 

AMS subject classifications. 


1. Introduction. Since its introduction by Bellman in the 50’s, dynamic pro¬ 
gramming has become a fundamental tool in the design of optimal control strategies 
for dynamical systems. It characterizes the value function of the corresponding op¬ 
timal control problem in terms of functional relations, the so-called Bellman and 
Hamilton-Jacobi-Bellman (henceforth HJB) equations. We begin by briefly recalling 
this setting in the context of infinite horizon optimal control. 

We make the following assumptions. We equip for n G IN with the Euclidean norm 
II • II 2 . Furthermore, let 

\f{x,u)-f{y,u)\<ujR\\x-y\\ 2 , 
\lix,u)-l{y,u)\<L0R\\x-y\\2, 


for x^y e with \\x — y \\2 < R and modulus ujr: [0, oo) ^ [0, oo) of polynomial 
growth satisfying lim^^o+ = 0 (cf. Ishii [H]), d,m G IN. Let the dynamics be 

given by 

O', / y(t) = 

\yio) = x 

for i > 0, where a; € and u G U = {u: M+ U measurable}, U C compact. 
We introduce the following cost functional J :V( 

pOO 

J{u)= / l{y{s),u{s))e~^^ds, A>0, 

Jo 

where y is the solution of ( |1.2[ ) depending on x and u. By the application of the 
dynamic programming principle, the value function 


v{x) = inf J{u) 

u^lA 
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is characterized as the viscosity solution [H Chapter 3] of the HJB equation 

Xv{x) + sup{—/(x, u) • Vv(x) — l(x, u)} = 0, X e 
ueu 

There exists an extensive literature concerning the construction of numerical schemes 
for static HJB equations. The spectrum of numerical techniques includes ordered 
upwind methods EHIE], high-order schemes m, domain decomposition techniques 
[7] and geometric approaches [6], among many others (we refer to [131 Chapter 5, 
p.l45] for a review of classical approximation methods). In this paper we follow a 
semi-Lagrangian approach m, which is broadly used to approximate HJB equations 
arising in optimal control problems, see,e.g., [13]. We illustrate the basic steps in the 
formulation of a semi-Lagrangian scheme for our model problem. 

The construction of a first-order semi-Lagrangian scheme begins by considering 
an Euler discretization of the system dynamics with time step h > 0 

( y-+^=y- + hfiy^,U^), 


for n G IN^, x G and controls G U. Then, the application of the dynamic 
programming principle on the discrete-time dynamics leads to the Bellman equation 

v{x) = min{(l — \h)v{x + hf{x^ u)) + hl{x^ u)}^ x G 

uEU 

To discretize this equation in space we introduce a bounded domain ft = [a,b]^ C 
a^b G M, where we define a regular quadrangular mesh with N nodes and mesh 
parameter k. We denote the set of nodes hy ft^ C ft. Let the discrete value function 
be defined in all grid points, V := . However, note that x + f{x^u) for 

X G ftis not necessarily a grid point, and therefore the value function has to be 
evaluated by interpolation which is chosen as a linear one here. The interpolant 
I[']{x) is defined on the basis of the dataset V. The resulting fully discrete scheme 
then reads 

(1.3) V{x) = min{(l — Xh)I[V]{x + hf{x, u)) + hl{x, rt)} =: G(E), x G ftk, 

u£U 

which can be solved by a fixed point iteration starting from an initial guess by 


An alternative to the fixed point approach is the use of Howard’s algorithm or iteration 
in the policy (control) space [5] [T], which is faster but also utilizes the parametric 
minimization of the Hamiltonian. 

Finally, once the value function is computed this allows to derive a feedback control 
for a given state x G by 

u{x) G argmin{(l — Xh)I[V]{x + hf{x, u)) + hl{x, i^)}. 
ueu 

A characteristic feature of the class of HJB equations arising in optimal control 
problem is its nonlinear Hamiltonian, 

sup{—/(x, u) • Vv(x) — l(x, 1 ^)} , 
ueu 
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which requires a parametric maximization (minimization) over the control set U. For 
its discrete analogue G{V), a common practice in the literature is to compute the min¬ 
imization by comparison, i.e., by evaluating the expression in a finite set of elements 
of U (see for instance [U [12] |T7| and references therein). In contrast to the com¬ 
parison approach, the contribution of this paper is to demonstrate that an accurate 


realization of the min-operation on the right hand side of (1.3) can have an impor¬ 
tant impact on the optimal controls that are determined on the basis of the dynamic 
programming principle. In this respect, the reader can take a preview to Figure [4^ 
where differences between optimal control fields obtained with different minimization 
routines can be appreciated. Previous works concerning the construction of minimiza¬ 
tion routines for this problem date back to [8], where Brent’s algorithm is proposed 
to solve high dimensional Hamilton-Jacobi -Bellman equations and to m, where the 
authors consider a fast semi-Lagrangian algorithm for front propagation problems. In 
this latter reference, the authors determine the minimizer of a specific Hamiltonian by 
means of an explicit formula. Moreover, for local optimization strategies in dynamic 
programming we refer to m for Brent’s algorithm and to [22] for a Bundle Newton 
method. 

In this article, a first-order primal-dual method (also known as Chambolle-Pock 
algorithm i) and a semismooth Newton method [ISl[2l] are proposed within the semi- 
Lagrangian scheme for the evaluation of the right hand side in (1.3). In contrast to 


the minimization by the comparison approach, the proposed algorithms leads to more 
accurate solutions for the same CPU time. Since we preserve the continuous nature 
of the control set, in some specific settings it is also possible to derive convergence 
results for our minimization strategies. Furthermore, it provides a solid framework 
to address challenging issues, such as nonsmooth optimal control problems with £i 
control penalizations in the cost functional. 

The paper is organized as follows. In Section we begin by recasting the dis¬ 
cretized Hamiltonian as a minimization problem explicitly depending on the control 
u. In Sectionwe introduce and adapt the Chambolle-Pock and semismooth Newton 
methods for the problem under consideration, and in Section we present numerical 
examples assessing the accuracy and performance of the proposed schemes. 

2. Explicitly control-dependent Hamiltonians. In this section we study the 
numerical evaluation of the discretized nonlinear Hamiltonian 


(2.1) min{(I — Xh)I[V]{x + hf{x, u)) -1- hl{x, u)} . 

uEU 

This is a non - standard minimization problem requiring the evaluation of the nonlin¬ 
ear mapping u I[V]{xhf{x,u)) which depends on the discrete dataset V and the 
system dynamics /(x, u). Therefore, a first step towards the construction of local min¬ 
imization strategies is to recast ( |2.I| ) by assuming specific structures for /[U],/(x, i^), 
and l{x,u)^ leading to explicit piecewise linear or quadratic optimization problems on 
u. This section is split between the treatment of the interpolation I[V]{x + h/(x, u)), 
and the evaluation of /(x, u). For the sake of simplicity we restrict the presentation to 
the two dimensional case d = 2, although the presented framework can be generalized 
to higher dimensions. 

2.1. The interpolation operator and local subdivisions of the control 
space. To analyze the interpolation operator we consider for every node x = (xi, X2) G 
flk the patch of four triangles defined by the neighboring nodes, cf. Fig. |2.I[ Assume 
that the arrival point x-\-hf{x, u) for x E u E U, is located in a triangle defined by 
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the points x = ^ ^ and x^ in with associated values Vi, V 2 , and V 3 as indicated 

in Fig. |2.1[ The linear interpolation formula then reads for a mesh point x G flk 

( 2 . 2 ) I[V]{x) = cxi + dx 2 + e 


with 



Fig. 2.1 : Arrival points and related control sets. 


Remark 2.1. Note that in the considered case of first order approximations the 
passage from the interpolation on one triangular sector to another one is continuous. 
For the first order interpolation schemes in a point x G one might also think 
of considering one interpolant for a macro-cell defined by a larger set of neighboring 
nodes. This could be computed by using 2^ suitable points. However, in practice 
this is similar to consider a grid of size 2 h, which leads to less accuracy, which is, in 
particular, in higher dimensions a crucial issue. Alternatively one could consider one 
interpolant for the macro-cell which is defined piecewise over all the triangulars. How¬ 
ever, this piecewise smooth interpolation is not convenient for numerical optimization. 
As a further alternative one can resort to higher order interpolation schemes. □ 

According to (113, the interpolation operator I[V]{x) is evaluated at the arrival 
point X -h hf{x,u). If this point is sufficiently close to x G O/., it is in general 
possible to know a priori in which of the four triangular sectors it will be located. 
Our goal is to establish a consistent procedure to locally divide the control space in 
sets that generate arrival point related to a single triangular sector. For instance, 
in the case of the eikonal equation there is a simple correspondence between the 
computed triangulars in control and state spaces, see Section [23) In the following we 
restrict ourselves to dynamics of the form 

(2.4) f{x,u) = g{x) + Bu, 

for B G i.e. nonlinear in the state, and linear in the control. We assume that 

h > 0 is chosen such that 
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Let the rows of B be denoted by {bi}f^i. Further, let xd denote a grid point in 
let X C {1,..., d} be an index set with complement and define a triangle Qx C 
associated io xd for the interpolation and X by 

Qx = {x = Xd x \ ||x||i <k,Xi>^ for i G X, < 0 for i G X^}. 


Then 


Ux = {u \ g{xD)i + biU >0 for i G X, g{xD)i + biu < 0 for i G X^} 

has the property that xa = xd + h{g{xD) + Bu) G Qx for u G Ux, and 0 < h < h. In 

) that ||x||i < k, Xi > 0 

for i G X and < 0 for i G X^. Note also that 

u= [jUi, 

Tew 


fact, X = Xa — Xd = h{g{xD) + Bu)^ and we obtain with (2.5 


where W is the set of all subsets of index sets in {1,..., d}. The associated interpola¬ 
tion operator on Qx is denoted by Jx- Note that evaluation of Ix[V]{x + f{x,u)) for 
X G Qx and u eUx leads to a linear dependence on u of the form 


(2.6) 


cx^i + dxU2 + ex , 


where the coefficients cj, dx and ex are uniquely determined by the vertices of Qx- 

2.2. Evaluation of the cost term. Having approximated the interpolation 
term of the Hamiltonian, what is left is to provide an approximation of the running 
cost l{x^u). If this term is defined in a pointwise manner, this imposes no additional 
difficulty. For instance, in some of the examples we shall utilize 

(2.7) K^^u) = \\x\\l + '^\\u\\l=x\+xl + '^{u\+ul) x G fife, ueU. 


This introduces a constant and a quadratic term in u to be added to the above 


presented expressions for the interpolant (2.6). For I given as in (2.7) we will consider 
in the scheme (m the explicit dependence on u. For a given node x G we have 


(2.8) 


[Vl = mm{Hx{x,V,u,I)}, 


where [V]x denotes the value of the discrete value function at node x. Further, for 
each set X we have 


(2.9) Hx{x,V,u,l) 


( 2 . 10 ) 


min {pix\y]{x + hf{x, u)) + hl{x, i^)} , 
ueUx 


mm 

ueUx 


cxUi + diU 2 + ex) + h (x^ + 2/^ + ^ (^ 

dx 2 , 2 


min {—Ml + —U2 + ciui + dxU2 + ex} , 

ueUx 2 2 



with 


/d = 1 — Ad , dx = dy , 6 x = dy , cx = Pcx , 
dx = pdi , ex = pex + h{x^ + y^). 
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Remark 2.2. For minimum time optimal control problems, the resulting semi- 
Lagrangian scheme for points x G Qk reads 

[V]x = mm{^I[V]{x + hf{x + u)) + 1 — , where (3 = , x G 

uEU 

see m, and the Hamiltonian takes the simplified form 

Hx{x,V,u,I) = min{cxtii + JxR 2 + ex}. 
ueUx 


□ 

To evaluate the right hand side in (2.8) numerically, we compute (2.9) for X G W by 
applying numerical optimization methods and then determine the minimizer on the 
macro-cell by comparison over the sets Ux- 

To determine the minimizer of the right hand side in (2.10) we define 


^ + cxui + dxU2 + ex- 

Then the minimization problem is given by 


( 2 . 11 ) 


min F{u), I e W. 

uEUx 


In particular, the following types of cost functionals are of interest for scalars 
a, 6 , s G M, where a, c > 0. Thereby we include also a bilinear term buiU 2 

to allow also bilinear interpolants over quadrangular cells. For the triangular inter¬ 
polation used throughout this paper, we take 6 = 0 . 

1 . Infinite horizon problem with quadratic control cost: 

(2.12) -F2(ri, U 2 ) = —+ buiU2 + —^2 3 - dui -h eu2 + t . 

The functional is smooth and convex if ac — 6 ^ >0. 

2. Minimum time problem: 


(2.13) 


Fmt{ui,U 2 ) = buiU 2 + dui -h eu 2 + r. 


The functional is non-convex if 6 7 ^ 0 but smooth. 

3. Minimum time/infinite horizon problem with ^i-control cost: 


(2.14) 


^mt/ih('^15'^2) = Ci\ui\ + buiU2 + c\u2\ + dui + eu2 + r. 


The functional is non-convex if 6 7 ^ 0 and non-smooth. 

4. Infinite horizon problem with ^ 2 - and ^i-control cost: 

(2.15) FiYi{ui,U 2 ) = + l\ui \ + buiU 2 + ^ul + s\u 2 \ + dui + eu 2 + r. 

The functional is convex if ac — 6 ^ >0 and non-smooth. 

Remark 2.3. Since U is compact all four functionals allow a minimum. If 


ac — 6 ^ > 0 it is unique for (|2.12[) and (|2.13[). For a discussion of optimal control 


problems of dynamical systems for functionals (2.14) and (2.15) we refer to [anH]. □ 
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2.3. Special case: eikonal dynamics. We consider the relation between the 
subdivision of the state and control space in the context of minimum time problems 
with eikonal dynamics. To set up the control problem we introduce a closed target 
T C , with int(T) 7 ^ 0, and smooth boundary. We consider the minimum time 
function 


(2.16) t{x^u) 


inf {t G M+ : ^(t, u) G T} if ^(t, u{t)) G T for some t, 
+OC otherwise, 


where u ^ hi and denotes the solution of ( 1 . 2 ) depending on the control. 

Furthermore we assume a small-time local controllability assumption as formulated 
in [13 p. 216]. The minimum time problem is defined as 


T{x) = inf t{x,u) 

u^U 

and the minimum time function can be characterized as the viscosity solution of 
sup {—/(x, u)^\/v{x)} = 1 in 7^ \ T, 

(2.17) 


v(x) =0 on ST, 


where 7Z are all points in the state space for which the time of arrival is finite. 
We consider the dynamics 


(2.18) f(x, w) = (^ “1 j , U = {uGR^: ||w ||2 < 1}, a: e 

which leads to an equation of eikonal type. For this problem we have a direct relation 
between the direction of the control u and the identification of the arrival point in the 
state space; for instance, all the arrival points in the first quadrant will correspond 
to controls belonging to the set U^i_ 9 } = G | 0 < , 0 < 1^2 }. It is possible to 

express the interpolation operator ( |1.3[ ) for given node x piecewise as 

(2.19) I[V]{x + hf{x,u)) = Ix[V]{x + hf{x,u)) 

with X + hf{x^u) in quadrant For every quadrant we obtain control-dependent 
formulas of the form 


( 2 . 20 ) 


Ix\V]{x + hf{x,u)) = cxui + dxU 2 + ex 


with coefficients cx^dx^ex G 


3. Algorithms for solving minimization problems of the form (2.11). We 

consider the following approaches: 

a) Minimization by comparison over a finite subset t/finite C Ux- 

b) First order primal-dual method: Chambolle-Pock algorithm, see [9]. 

c) Second-order method: Two different types of semismooth Newton methods 
depending on the smoothness of the cost functional. 


d) If the functional is of type (2.13), the controls are located on the sphere or 
in the origin and can be found by a classical Newton method with a suitable 
chosen initialization over the parameterized sphere. 
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The minimization by comparison approach a) is broadly used in the literature. It 
consists in choosing a finite subset t/finite U where the cost function is evaluated and 
the minimum is selected among the corresponding values. Such a procedure induces a 
different optimization paradigm, in the sense that the continuous nature of the control 
space is replaced by a discrete approximation. If a parameter-dependent discretization 
of the control set is considered, where errors with respect to the continuous set can 
be estimated, the discretization has to be fine enough such that the error is negligible 
compared to the errors introduced in the different discretization steps of the overall 
scheme. Furthermore, accurate discretization of the control space are far from trivial 
even in simple cases, as in the three-dimensional eikonal dynamics with U = G : 
ll^lb < 1}, where a spherical coordinate-based discretization of the control introduces 
a concentration of points around the poles. In contrast to this approach, we propose 


several numerical algorithms for the treatment of problems of the form (2.11), which 


preserve the continuous nature of the optimization problem at a similar computational 
cost. 


3.1. The smooth case: Cost functionals of type ( |2.12 ) and ( 2.13[ ). In this 
section we consider control constraints of the type U = {u ^ • ||'?^ 2 < !}• To 

solve the optimization problem for smooth cost functionals of type (2.12) and ( |2.13 ) 


we present a first-order primal-dual algorithm and a semismooth Newton method. 


3.1.1. Primal-dual algorithm. 

reformulated as 


Here we assume that problem (2.11) can be 


(3.1) min F{u) + Ik{u), 

where F is smooth, and convex and Ik{u) is the indicator function oi K = Ux- This 
fits in the setting presented in [9], where a primal-dual algorithm (also known as 
Chambolle-Pock algorithm) is formulated, when using the same F, and by choosing 
the mappings denoted by K and G in the reference by id and Ik- We recall the 
algorithm in the following. 

Algorithm 1: Chambolle-Pock algorithm 

Data: Choose n = 0, r > 0, cr > 0, > 0, 'd G [0,1], {uo,yo) ^ x 1^^, 

uo = Uq. 

repeat 

Compute Xn+i = (^n+i,f^n+i) by 

Vn+l = (^ + CrdF*)“^(^n + CrUn)^ 

Un^i = (^ + rdG)~^{Un - r^n+l), 

^ "^n+l ~ "^n+l T "^n)* 

Set n = n 1- 

until \\xyi 1 IliRd V'-! 

From identities from convex analysis, see m, we have 

(3.2) (J + TdGy{y) = argminji ^^^—— + 

uGM™ I ^ J 

u-(/ + raF)-i(w)=r(/+iaF*) (“) , 


(3.3) 
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where Pk denotes the projection on K. From (3.3) we have with cr = ^ 

-1 


and hence 


(/ + adF*) ^{au) = au — a [ I —dF ) (i^), 


(j + (TaF*)-i(u) = u-(T(/+-aFj (-u). 


a j \G 

In the case X = {1, 2} the projection Pk is given by 

max(0,p) 




K • 


Pk{p) = 


max(l, ||max( 0 ,p)|| 2 )’ 


This leads to Algorithm 


Algorithm 2: Chambolle-Pock algorithm (variant) 


Data: Choose n = 0, r > 0, cr > 0, > 0, 'd G [0,1], {uo,yo) e x 

Uq = Uq. 

repeat 

Compute Xn+i = (2/„+i, m„+i, itn+i) by 


Vn+l — Vn T ^ 




'^n+l — PRij^n '^yn-\-l)s 

< — ^n+1 T ^n)* 

Set n = nPl. 
until ^Xji X 72—1 lli^cixR"*'xR"*' ^ V'! 


If (/ + ^dP) ^ is linear, the first step can be reformulated as 

-1 


Vn+I = (^I+ OF + «„) . 


Remark 3.1. Since the cost functional has only quadratic and linear terms it 
is more convenient to use a linear interpolant rather than a bilinear one in the semi- 
Lagrangian scheme.This does not produce any bilinear terms. If F has a bilinear term 
it is very challenging to compute (/ + ^dF)~^ in higher dimensions, since then dP 
depends nonlinearly on u. □ 

3.2. Semismooth Newton method. The Chambolle-Pock algorithm is a first 
order algorithm which requires convexity of the functional. We refer to [23] for an 
extension to a class of non-convex problems. As a second algorithm that we propose 
we turn to the semismooth Newton method which does not rely on global convexity. 
We recall some main aspects of semismooth functions, cf. [21]. 

Definition 3.2 (Semismoothness). Let V C be nonempty and open, /i G IN. 
Then funetion f: V ^ is semismooth at x ^ V if it is Lipsehitz eontinuous near 
X and if the following limit exists for all s , z/ G IN.* 

lim Md. 
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Furthermore, there holds the following chain rule. 

Lemma 3.3 (Chain rule). Let V C and W G ML be nonempty open sets, 

g: V ^ W he semismooth at x E V, and h: W ^ 77 G IN, be semismooth at g{x) 

with g{V) C W. Then the eomposite map f := h o g: V ML is semismooth at x. 
Moreover, 


fix,-) = h'{g{x),g\x,-)). 


To set up a semismooth Newton algorithm for (2.11) we proceed as follows. For 


simplicity, we focus on the case X = {1,2}. The first-order optimality condition for 
cost functionals of type ( 2.12[ ) and ( 2.13| ) can be formulated as 


(3.5) 


u = Pfc( 7 x — 'dVF(u)) V'd > 0 


for u G F = F 2 or F = Fmt with projection Pk as in (3.4). We introduce 


p = y — 'dVF{y) and choose id such that idVF{u) is of the same scale as y and rewrite 
(|3.5|) as 


u - Pk{p) = 0, 
u — id\/F{u) — p = 0. 

Setting (3 = max(l, || max(0,p)||2) we define 

( pu — max(0,p) 
u — idVF{u) — p 
^-max(l, ||max( 0 ,p)|| 2 )y 

Now the optimality condition can be formulated as 

E{z) = 0 

with z = {u,p,/3). Then we set up a semismooth Newton method as presented in 
Algorithm 

Algorithm 3: Semismooth Newton algorithm 

Data: Choose n = 0, initialization zq e Ux x x M, 77 > 0. 

repeat 

Solve the Newton equation for 5z £ Ux x x M given by 


(3.6) 


Pnl 

I - dV^F{Un) 


D 




0 




Xpn >0 

-I 

-pId,, 


yn\ 

0 

ly 


5 z — F ( Zji ) 


(with matrix = diag(xp^>o) and rUn = || max( 0 ,Pn)|| 2 )- 

Update 


(3.7) 

Set n = n 1. 

until \\Zn - Zn-l\\l 


Zn-f'i- — ^Ti T 


< m 
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We address solvability of (3.6). 

Lemma 3.4. Set M = I — ^V‘^F{un) and = max(0,g') for q G 1 
Newton matrix is regular if one of the following two conditions is satisfied: 
1. rrin < 1: /?„ ^ 0 , f regular. 


The 


2. rUn > I: /?„ ^ 0, 


(3.8) 


I — —M regular 

Pn 


and |-^(p+)r(J-^M+)-i (-M+A+ 4 ))+ + i| ^q, 


where the 


components of p = {p'^,p~) are ordered with respect to active and inactive 
sets and M+ denotes the submatrix corresponding to . 

Remark 3.5. For ~ 1 condition (3.8) is closely related to the regularity of 
the Hessian of F. Moreover, if we set K = ||V^F(i 4 )|| and assume that Pn 

are bounded away from 1 (i.e. /3n > /^ > 1) then the choice'd < implies (3.8). □ 

Proof of Lemma \3.4\ The two cases are considered separately. For mn < 1, the 
Newton-Matrix is given by 


^Pnd ^Xpn> 
DE{Zn) = I M -I 


0 


0 


The 


regularity of the matrix follows from the regularity of ^^d 

(;^^^Xpn>o + • 

For uin P 1 the Newton-matrix is given by 

^Pnl ^Xpn>0 

DE{Zn) = { M -I 0 
0 1 


hence 


and leads to the Newton equations 


(3.9) 

(3.10) 

(3.11) 


I3n5u - (5p+ + 6l3Un = Zi, 
MSu — 6p = Z 2 ^ 

— —ip'^fSp++5p = Z3 


for suitable zi,Z 2 G 2(3 G M independent of {6u, 6p, 6P). From the first and second 
equations we obtain 

- Tm((5p+ - 5Pun) +5p = Z 4 

Pn 

for 2(4 G Without loss of generality we assume that the components of 6p = 

{6p'^,6p~) are ordered with respect to active and inactive sets. Then we have 

1 11 ^ 

(/ - —M+)5p+ + SI3M+^ = Z 4 

Pn Pn 


(3.12) 
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and therefore 


(3.13) 



With (3.11) we obtain 


(3.14) 



Am+)“^ (-m+^ 

Pn V Pn 



+ 1 


5P = P. 


□ 


Theorem 3.6. Let u be a striet loeal minimizer of (2.11). Under the assump 


tions of Lemma \3.4\ the semismooth Newton method eonverges loeally superlinearly or 
terminates after a finite number of steps at u. 

Proof The operator E is semismooth and is as a composition of Lipschitz contin¬ 
uous functions again Lipschitz continuous. Furthermore from Lemma [3.4| we obtain 
the boundedness of the inverse derivative in a neighbourhood of u. 

Consequently, the assertion follows from [HI p. 29] and [151 P- 220, Thm. 8.3]. □ 


3.3. Approach for cost functionals of type ( |2.13 ). To treat the problem 
( 2 . 11 ) for cost functionals of type (2.13) we present an alternative approach. We 


show that all possible minimizers are located on the sphere or in the origin; cf. also 
the discussion in m where a linear functional is considered. To minimize over the 
sphere we parameterize the sphere by polar coordinates and consider the restriction 
of the functional on the sphere. There the minimizer can be found by applying a 
classical Newton method. 

Lemma 3.7. For eontrol problems (2.11) with eost funetional (2.13) forr = 0 all 


minimizers are loeated either on the sphere or in the origin. 


Proof For simplicity we consider the case ( 2 . 11 ) for d = 2 and X = { 1 , 2 }. Let 
F{u) = dui + buiU 2 + eu 2 , rt G d, 6, e G M. Then we have 


\/F{u) = 


d + bu2 
e + bui 


V‘^F{u) = 


This implies, we have for 6 7 ^ 0 a saddle point with eigenvalues Tl in 

e d 

u, = --, u, = --. 

For 5 = 0 the minimum is reached in u = (0,0)^. Consequently, all minimizers are 
on the boundary of However, since the restriction of the bilinear functional 

to the axes is linear the assertion follows immediately. □ 

In two dimensions the functional F restricted to the sphere is given by 


~ 5 TT 

F{ip) = acos((p) + - sin(2(/:?) + csin((/:?), 0 < < —. 


3.4. The non-smooth case: Cost functionals of type (2.14) and (2.15). 


In this section we consider two types of control constraints, box constraints as well as 
Euclidean constraints. Finally, we also present a splitting approach. 
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3.4.1. Semismooth Newton in case of Euclidean norm constraints. In 

this section we restrict the consideration to equations of eikonal type and present 
the approach for a two dimensional problem. Therefore the control set is given by 
t/ = {i4GM^:||i4||2<l}. Let us consider the equivalent problem formulation 


(3.15) 


min Eiu) + hiu), 

neM2 ^ ^ ^ 


where E’: ^ M is smooth and 


(3.16) 


h{u) 


a{\ui\ + \u 2 \), Hue Lx, 
oo, H u ^ Ux- 


The optimality condition for the nonsmooth problem (3.15) is given by 


0 G \/E{u) + dh{u) , 


where dh denotes the sub differential of the convex function h. Setting q = —VE{u) 
the optimality condition can be equivalently written as 


(3.17) 


{q = -WE{u), 
^ g e dh{u). 

With the convex conjugate we obtain equivalently 


(3.18) 

with 


q = -VE{u), 
u G dh*{q), 


Note that on Ux 


h*{q) = sup {v ' q — h{v)). 
veUx 


v-q- h{v) = {qi - ai)vi + {q2 - a2)v2 

with 

(3.19) ai = sgn(Li)(a, 0^2 = sgn(L 2 )Q^ 
and 

(3.20) 

_ J 1 , ifl= { 1 , 2 } or 1 = {!}, r _ / 1, if J= {1,2} or / = {2}, 
^ \ -1, ifl= {2} orl = 0, ’ ^ \ -1, if/ = {!} or2: = 0. 

For fixed q and 1 ; G Lx we define 

Kv) = {<li - «i)vi + {q2 - a2)v2. 

The sup of this linear function is attained on dUx- We distinguish four cases 
sgn(gi - ai) 7 ^ sgn(Li) A sgn(g2 - a2) 7 ^ sgn(i2) : hl^ L^{q) = 0, 

sgn(gi - ai) = sgn(Li) A sgn (92 - ^2) 7 ^ sgn(L2) : hl^ j^^{q) = qi - a, 

sgn(gi - ai) 7^ sgn(Li) A sgn(g2 - ^2) = sgn(L2) : hl^ j^^{q) = q2 - a, 

sgn(gi - ai) = sgn(Li) A sgn (92 - a2) = sgn(i2) : 

^LulM) = sup(gi -Q!i)cos(/?+ {q 2 -a 2 ) sirup, 

(peA 
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with 


(3.21) 


A = 


^ 0 < (^ < f ifX= {1,2}, 
\ <'K \iX — {2}, 

TT < ip < ^ ifX=0, 

^ ^ < ip <2 ti if X = {!}. 


In the following for simplicity we only consider the problem for X = {1, 2}. From 
the first order condition we obtain the following relation between q and ip 


(3.22) 

Summarizing we have 


tan ip = - > 0. 

Qi-a 


(3.23) 


h*{q) = 


(gi - (a) cos ip + {q 2 - <a) sin ip, 

if qi,q2 > O', 

<72 - 

if qi < a, q 2 > a 

0, 

if qi,q 2 < o, 

gi - a. 

ii q\> a,q 2 < o, 


with ip given by (3.22) and for the generalized derivative 


(3.24) 


where 


dh*{q) = < 


{wi{q),W 2 {q)P, if qi, q 2 > a, 

(0,1)^, if Qi <a, q 2 > a, 

(0,0)^, ifqi < (a, q2 < (a, 

, (1,0)'^, if qi> a, q 2 <a 


{wi{q),W 2 {q)) = V{{qi - a) cos(^ + (^2 - a) sin(/?). 

To set up a semismooth Newton method we introduce a piecewise affine, contin¬ 
uous regularization of the gradient in a 5-tube around the discontinuity and define 


(3.25) 


' {wi{q),W2{q)P, 
( 0 , 1 )^, 


{dh*)e{q) = < 


q2-a 


( 0 , 0 )^, 

qi-a 


( 1 , 0 ) 
r 
s 


,0 


if gi, q2> 

if gi < (a, q 2 > a-\- 5, 
if gi < (a, a < q 2 < a e, 
if qi < a, q2 < a, 
if a < qi < a s, q 2 < a, 

if gi > (a -h 5, q2 < (a, 


-{wi{qp),W 2 {qp))^, if a <qi, a < g2, and ||g - ((a,(a )^||2 < 5, 


where 


^ = Ik “ (<a, <a )||2 5 ^ = tan ^l ——^ ) , and = 5(cos(^), sin(^)) + (a, (a). 

' qi-a ’ 
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The semismooth Newton method is presented in Algorithm 


Algorithm 4: Semismooth Newton algorithm 

Data: Choose n = 0, initialization zq = (qo,uo) G x Ux, rj > 0. 

repeat 

Solve the Newton equation for G x Uj given by 


(3.26) 


f I V^E(Un)\ c _ f <ln + VE(Un) \ 

\-D{dh*Uqn) I ) \-{dh*)e{qn) + Un) 


where D{’) denotes the Newton derivative. 
Update 


(3.27) 




Set n = n 1. 

until \\Zn - 2;n-l||R2xR2 < T]; 


We can conclude the following theorem. 

Theorem 3.8. Let u be a striet loeal minimum of (3.15). If 


(3.28) 


I-V^E{u„)D{dh*Uqr,) 


is regular for all k, the semismooth Newton eonverges loeally superlinearly. 
Proof We introduce the piecewise affine function 


G: X ^ X G{q, u) = [q^ + VE{u)^ — {dh^)s{q) + Un)- 


G is semismooth and together with condition (|3.28| ) we get directly the regularity of 
the Hessian in the Newton equation given in (3.26). From both conditions we derive 
locally superlinear convergence as in Theorem |3.6| □ 


3.4.2. Semismooth Newton in case of box constraints. Now, let the con¬ 
trol set be given hy U = G BN : Ua < u < ui,} with ^ 0 ^ '^6 ^ B‘^. As 

in the previous section, we consider a cost functional of the type 


(3.29) 


min E(u) h(u), 

neM2 ^ ^ ^ 


where E is smooth and 
(3.30) h{u) 


a(|wi| + |m2|), iiueUx, 
oo, if u ^ Ux- 


Proceeding as in the previous section we distinguish four cases 


sgn(gi - ai) ^ sgn(Li) 
sgn(gi - ai) = sgn(Li) 
sgn(gi - ai) ^ sgn(Li) 
sgn(gi -ai) = sgn(Li) 


A sgn(g2 - ^ 2 ) sgn(L2) 
A sgn(g2 -a2) sgn(i2) 
A sgn(92 - ^ 2 ) = sgn(L2) 
A sgn(g2 - ct2) = sgn(L2) 


— 0 ; 

= kl - Oii\Ua, 

= 1^2 - a2|M6, 
= ki -ai\ua 
+ \q2 - a2\Ub , 
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with ^ 1,^2 as in (3.19) and as in (3.20). Further, 


there holds 


(3.31) 


Vh*{q) 


{sgn{qi) Ua,sgn{q 2 ) UbX, 

if 

Qi 

> 

a. 

q2 

> 

a. 

{0,sgn{q2) Ub)'^, 

if 

qi 

< 

(a. 

q2 

> 

(a. 


if 

qi 

< 

a. 

q2 

< 

a. 

{sgn{qi) Ua,0X, 

if 

qi 

> 

a. 

q2 

< 

a , 

ingly for the other cases. 

To 

set up 

a semismooth Newton 


method we regularize the first component of the gradient with ramps of width s at 
\qi\ = a and the second component with ramps at \q 2 \ = a. The overall algorithm 


has the same form as the semismooth Newton method in (3.26). 


By a similar consideration as in the previous section Theorem 3.8 holds also for 
the semismooth Newton method described above. 


3.4.3. A splitting approach. To treat the cost functionals (2.14) and (2.15) 
with -terms we can reformulate the problem without non-differentiable terms by 
doubling the number of variables. This is illustrated for the two dimensional case 
with two dimensional control {u^v) G Ux- For a scalar z G M we define = 
max(0, z), Z- = min(0, z) (in particular =0 and \z\ = Z-^ — Z-). Then problem 


(2.11) with cost functional (2.15) is given by 


min au^ + au_ + cvt + cv_ + bu-^V-^ — bu-V-^ — bu-^V- + bu-V- 
+ — U-) + s{v-^ — V-) + d{uj^ + U-) + e(i;+ + V-) + r, s.t. 

Qi + h ("“+ “ ““i > 0, 

— V- J 

= 0 , = 0 , 

>0, 0, U- >0, V- > 0, 

+ ll^lll < 1> u,veUi 


with g chosen as in (2.4) and bi as the columns of B. This problem formulation can 


be simplified, by distinguishing four cases depending on X; for example, for F^{i, 2 } we 
obtain the equivalent problem 


(3.32) 


min au\ + cv‘^ + bu-^V-^ + (/ + d)u-^ + (s + e)v-^ + r, s.t. 

+ 1 = 1 , 2 , 

U-\- >0, V-^ > 0, 

v‘^ < 1. 


Problem (3.32) can be solved by using algorithms of Section 3.1 


4. Numerical examples. In this section we present a set of numerical tests 
aiming at studying the performance and accuracy of the algorithms presented in the 
previous sections. We begin by assessing the performance of the proposed numerical 
optimization routines as a separate building block. 
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4.1. Preliminary tests. We consider a generic two-dimensional minimization 
problem of the form 

(4.1) min ^\\u\\l + L ■ u + ^\\u\\i 

u£U Z 

subject to Euclidean norm or box constraints, i.e. U = {u G I \\u \\2 < 1} or 


constraints. For every setting we can observe that the minimization by the comparison 
algorithm (see Section requires very fine discretizations of the control variable 
(and thus, higher CPU time) to achieve similar error levels as the optimization-based 
counterpart. 


4.1 


U = {u |0<14i< 1,0<142<1}, respectively. Results presented in Tables 
|4.2| and |4.3| show the different performance scenarios found under different costs and 


Algorithm 

Tolerance 

Iterations 

CPU time 

• 2 error 

Chambolle-Pock 

lE-4 

9 

7.3E-5 [s] 

4.31E-5 

Semismooth Newton 

lE-4 

5 

1.3E-2 [s] 

7.74E-9 

Comparison (1E4 evaluations) 

- 

- 

l.lE-3 [s] 

1.6E-2 

Comparison {2E3 evaluations) 

- 

- 

6.6E-4 [s] 

4.1E-2 

Table 4.1: Performance tests for 

an ^2-cost (i.e. 

7 = 0 ) with Euclidean norm constraint. 


Algorithm 

Tolerance 

Iterations 

CPU time 

• 2 error 

Semismooth Newton 

lE-4 

101 

4.53E-3 [s] 

1.51E-3 

Comparison (4E4 evaluations) 

- 

- 

5.18E-3 [s] 

5.77E-3 

Comparison (2E5 evaluations) 

- 

- 

2.08E-2 [s] 

2.80E-3 


Table 4.2: Combined £2- and -£i-cost (7 = 0.1) with box constraint. 


Algorithm 

Tolerance 

Iterations 

CPU time 

• 2 error 

Semismooth Newton 

lE-4 

58 

1.47E-2 [s] 

1.23E-3 

Comparison {2E3 evaluations) 

- 

- 

8.56E-4 [s] 

4.18E-2 

Comparison (1E5 evaluations) 

- 

- 

1.20E-2 [s] 

1.47E-2 


Table 4.3: Combined £2- and ^i-cost (7 = 0.1) with Euclidean norm constraint. 


4.2. Interplay with the semi-Lagrangian scheme. Having embedded the 
minimization routines within a semi-Lagrangian scheme, we show in Fig. |4.1| the 
evolution of the average iteration count per gridpoint (at every fixed point iteration 
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of the semi-Lagrangrian scheme), for both a minimum time and an infinite horizon 
optimal control problem subject to eikonal dynamics. Note the difference in the 
evolution of the subiteration count depending on the considered control problem. In 
the minimum time problem nodes that have not received information propagating 
from the optimality front are still minimized with irrelevant information until they 
are reached by the optimality front, whereas in the infinite horizon case “correct” 
information is available for every node from the first iteration due to the presence 
of a running cost. In this latter case, the impact of the available information from 
the previous semi-Lagrangian iteration is similar to a warm start of the optimization 
routine. 


Minimum time problem Infinite horizon problem 



Fig. 4.1 : Subiteration count for 2D control problems with eikonal dynamics 


4.3. Infinite horizon problems with ^2 running cost. We present three dif¬ 
ferent numerical tests for infinite horizon optimal control problem with cost functional 
and running cost given by 

J{u,x)= / l{x{s),u{s))e~^^ds , and = K Ikll2 + ll'^IlL 72 > 0, A > 0 . 

Jo 2 2 

As common setting, the fixed point iteration is solved until 


||yn+l_yn||<ifc2^ „ = 1,2,3,..., 


where k stands for the space discretization parameter and the stopping tolerance for 
the inner optimization routine is set to 10“^. Further common parameters are A = O.I 


and 72 = 2, and specific settings for every problem can be found in Table 4.4 


Test 

n 

u 

h 

Test I 

[-1,1]^ 

IH|2<i 


Test 2 

[-1,1]^ 

IH|2<i 

\k 

h 

Test 3 

[-1,1]^ 

ll'ulloo < 0.3 


Table 4.4: Parameters for Tests 1-3. 
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Test 1: 2D eikonal dynamics. In this test we consider eikonal dynamics of 
the form 


(4.2) f{x,u) = {Ui,U2f, ||(Wl,W2)||2 < 1 • 


We study the accuracy and performance of the semi-Lagrangian scheme with dif¬ 
ferent minimization routines: discretization of the control set and minimization by 
comparison, a semismooth Newton method given in Algorithmic and the approach 
by Chambolle and Pock given in Algorithmic In order to make a fair study of the 
different routines, we choose the discrete set of controls for the comparison algorithm 
such that the CPU time for the semismooth Newton method and the Chambolle-Pock 


algorithm is almost the same. Table 4.5 shows errors in both the value function and 
in the optimal control between the exact solutions v and u and their numerical ap¬ 
proximations Vh and Uh for the different schemes. Errors are computed with respect 
to the exact solution of the Hamilton-Jacobi-Bellman equation, which is not readily 
available in the literature. Since it is useful for numerical investigations, it is provided 
in Appendix B. The results show that we have similar CPU time of the approaches 
and independently of the meshsize, the optimization-based schemes yield more accu¬ 
rate approximations of the value function and the associated optimal control than the 
approach based on comparison. Further results are shown in Fig. |4.2[ where we also 
consider nonhomogeneous eikonal dynamics 


(4.3) 


f{x,u) = il + Xx2>0.5ix))iUi,U2)'^ , ||(Wi,M 2)||2 < 1 


with Xa: 2 >o. 5 (^) corresponds to the indicator function of the set {x = (xi,X 2 ) \x 2 > 
0.5}. The figure shows that both approaches, the SL-scheme with a semismooth inner 
optimization block and the one with a comparison-based routine, lead to very similar 
value functions. This is clearly not the case for the optimal control fields depicted 
in rows 2 and 3 of Fig. |4.2[ Even by a post-processing step it would be difficult to 
obtain the results in the third row from those in the second row. 


Algorithm 

k = 

0.05, (40^ 

DoF) 

k = 

0.025, (80^ 

DoF) 

CPU time 

l|v-T4||i 


CPU time 

Ik-Vfelli 

\\u-Uh\\i 

Comparison 

63.52 [s] 

3.12E-2 

3.84E-2 

5.76E2 [s] 

1.96E-2 

1.74E-2 

Semismooth Newton 

77.25 [s] 

2.62E-2 

1.61E-2 

7.27E2 [sj 

1.36E-2 

7.21E-3 

Chambolle-Pock 

63.05 [sj 

2.60E-2 

1.42E-2 

5.77E2 [sj 

1.36E-2 

6.83E-3 


Table 4.5: Infinite horizon control of 2D eikonal dynamics. CPU time and errors for a 
semi-Lagrangian scheme with different minimization routines. The comparison algorithm 
was run with a discrete set of 1280 control points in every node. 


Test 2: 3D eikonal dynamics. In order to illustrate that our approach can 
be extended to higher dimensions, we consider the infinite horizon optimal control of 
three dimensional eikonal dynamics 


f{x,u) = {ui,U 2 ,U 3 )'^ , ||(Ml,1t2,1t3)||2 < 1 . 
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Value fmictioii 


Value funetiou 




0|)iiiual cuulrul Uj via si?iiii>iiuuotli algurtllmi 


1 

0.5 

0 

- 0.5 

-1 

-1 


Optimal control ui via semismooth algorithm 


Fig. 4.2: Infinite horizon control with 2D eikonal dynamics. Left: continuous dynamics. 
Right: discontinuous dynamics. 


Errors and CPU time are shown in Table 14.61 Since the values for the semismooth 
Newton and the Chambolle-Pock algorithms are similar, we only include the data 
for this latter one. With the Chambolle-Pock algorithm we obtain more accurate 
solutions of the control field for a similar amount of CPU time as the comparison 
approach. 


Test 3: Triple integrator with two control variables. In this test the 
dynamics are given by a triple integrator with two control variables subject to box 
constraints 


f{x,u) = {x2,xs + '^ 1 , 1 ^ 2 )^, < Ci, \U 2 \ a,6 G R. 
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k = 

0.1, (203 DoF) 

k = 

0.05, (403 DoF) 

Algorithm 

CPU time 

h-Vh\\i 

\\u-Uh\\i 

CPU time 

\\v-Vh\\i 

\W-Uh\\i 

Comparison 

93.76 [s] 

2.49E-2 

2.46E-2 

2.89E2 [s] 

1.02E-2 

1.89E-2 

Chambolle-Pock 

97.25 [s] 

9.92E-3 

2.07E-2 

2.01E2 [s] 

5.66E-3 

1.22E-2 


Table 4.6: Infinite horizon control of 3D eikonal dynamics. CPU time and errors for a 
semi-Lagrangian scheme with different minimization routines. The comparison algorithm 
was run with a discrete set of 5120 control points. 


The purpose of this example is to stress that the minimization strategy that we have 
introduced can be also applied to non-eikonal dynamics, where the correspondence 
between octants in the state space and the control field is not trivial. For the sake of 
completeness the control space decomposition procedure for this example can be found 
in Appendix A. In this particular case, every octant will have associated a different 
rectangular sector in the control space for its arrival points. Results for the value 
function, optimal controls and trajectories are shown in Fig. | 4 . 3 [ In the second row 
of this figure, we observe distinct differences between the optimal controls obtained 
from the Chambolle-Pock and comparison-based algorithms. These differences in the 
control lead to different approximated steady states as it can be seen in the left of the 
second row of Fig. | 4 . 3 [ To highlight also the effect of closed-loop control we carried 
out an experiment with additive structural and output noise. In the third row of 
Fig. | 4 . 3 [ the resulting states from open and closed-loop control can be compared. 

4 . 4 . Infinite horizon problems with combined ^2 and ^i-costs. We present 
two numerical tests for infinite horizon optimal control problems with cost function 
and running cost given by 


fO 

J {u, x) = 

Jo 


l{x{s),u{s))e 


l{u, x) = 



72 

2 


NlU7ilklli, 


72 > 0 , 7i > 0 , A > 0 . 


The stopping rule for the fixed point iteration is defined in the same way as in the 
previous set of examples. All the optimization-based routines have been solved with 
a semismooth Newton method as the inner block, with a regularization parameter 
5 = le — 3 . Further common parameters are A = 0.1 and 72 = 2, and specific settings 
for every problem can be found in Table 


Test 

0 

U 

h 

k 

Test 4 

[-1,1]^ 

11^112 < 1 


0.025 

Test 5 

[ 0 , 27 r ]3 

||w||oo < 0.3 

0 . 2 k 

0.2 


Table 4.7: Parameters for Tests 4-5. 
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i«osurl<ice V — 1 


V — 2 


Iso.siiiicice V — 



-1 "-1 -1 -1 



-1 "-1 


Effect on the state trajectories of the minimization routine Effect on the input signal of the minimization routine 





Fig. 4.3: Triple integrator with two controls, numerical results with k — 0.05. Top: different 
isosurfaces. Middle: differences in the optimal control lead to different trajectories. Bottom: 
the feedback approach leads to robustness with respect to noise in the dynamics. 


Test 4: 2D eikonal dynamics. This test considers the same two dimensional 
dynamics as in Test 1, the only difference being the inclusion of an ^i-term in the 
cost functional. For the inner optimization block we apply the semismooth Newton 
method presented in Algorithm Results are shown in Fig. |4.4[ where differences 
in the shape of the value function can be observed as the 71 weight increases. In the 
second row of Fig. |4.4[ the effect of sparsity on the first control component can be 
seen from the fact that it is identically zero in a band around the origin. Moreover, 
the width of the band increases with 71 . As for the value function, introducing the 
^i-term breaks its radially symmetric structure, as it is shown in the first row of 
Fig. |4.4[ 

Test 5: 3D car model. In our last test, the dynamics are given by a nonlinear 
3D car model with two control variables: 


f{x,u) = {ui cos(x 3 ),?xi sm{x 3 ),U 2 ) 


with ui G [— cJi and U2 G [—002^ 002]^ > 0, C (;2 > 0. For this problem we imple¬ 

ment the semismooth Newton algorithm with box constraints introduced in Section 
3.4.2 Results are shown in Fig. |4.5[ The first row shows same isovalues for different 
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Fig. 4.4: Sparse control of eikonal dynamics. Top: inclusion of an ^i-cost breaks the radially 
symmetric structure of the solution. Bottom: the sparsity of the optimal control translates 
into a zero band around of the origin, depending on the weight 71 . 


costs. The addition of an -^i-term shrinks the region of a given isovalue. The second 
and third rows depict the effect of sparsity in both control variables, creating regions 
of zero action. A direct consequence of this can be seen at the bottom of Fig. |4.5[ 
where the inclusion of the additional ^i-cost generates optimal trajectories which do 
not reach the origin due to the interplay between the control cost and the discount 
factor of the optimal control problem. 

Concluding remarks. We have presented a numerical approach for the solution 
of HJB equations based on a semi-Lagrangian discretization and the use of different 
local minimization strategies for the approximation of the corresponding numerical 
Hamiltonian. The numerical results show a more accurate resolution of the optimal 
control field at a similar computational cost as the currently used schemes. Further¬ 
more, the proposed approach can be also adapted to treat non-differential costs such 
as ^i-penalizations on the control. Since the numerical approximation of the Hamilto¬ 
nian constitutes one building block within the construction of approximation schemes 
for HJB and related equations, the idea of using local minimization techniques can be 
conveniently recast in similar problems, such as front propagation problems and dif¬ 
ferential games, and in related approximation techniques, like fast marching schemes, 
policy iteration methods and high-order approximations. 

Appendix A: Decomposition of the control space for Test 3. We make 
an explicit presentation of the decomposition of the control space of Test 3. The 
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dynamics is given by 

f{x,u) = {X 2 ,X 3 , \ui\ < a, |tX 2 | <6, a,b eR. 

For a given departure point x^ = (xf, x^^x^) G R^ and a sufficiently small we want 
to identify a relation between subsets Ux of the control space U = [—a, a] x [-6, b] 
and the location of the corresponding arrival points 

x"^ + h{x2,X3 + Ui,U2) , {ui,U2) G Ux- 

In the three-dimensional case, the control set U can be decomposed into at most eight 
disjoint subsets, one per octant with 

u= [jUi 

XCW 

for W = {1,2,3}. Since every octant defines a unique linear interpolant, we solve 
eight different minimization problems, and then compute the global nodal minimizer 
by comparison. For instance, let us consider the sector 

Q{i,2,3} = {(a:i,a;2,a;3) \x^ < xi, < X2, xf < X3, {xi-xf)+{x2-X2)-\-{x3-xf) < k}, 

where the evaluation of the arrival point is defined by the linear interpolant /{ 1 , 2 , 3 } 
depending on the points x^ + (/c,0,0), +(0,/c, 0), and + (0, 0,/c). The subset 

^{ 1 , 2 , 3 } related to this interpolant reduces to 


= {{ui,U2) e c/1 rcg + Ml > 0, 1t2 > 0} . 


Note that for this definition to be meaningful, we need to assume that X 2 > 0, 
otherwise 2 , 3 } = 0- Also note that the condition (xi—xf) + (x 2 —X 2 ) + (x 3 —a; 3 ) < k 
is omitted by assuming a sufficiently small h ensuring it. The next subset ^^{ 2 , 3 } relates 
to the sector 

Q{ 2 , 3 } = {(X 1 ,X 2 ,X 3 ) \xf > Xi, X 2 < X 2 , xf < X 3 , -(Xi-X^) + (X 2 -X 2 ) + (X 3 -X 3 ) < k} , 
and is given by 


C^{2,3} = {{ui,U 2 ) eU\xf + Ui>0,U2>0}. 

Note that t/{i, 2 , 3 } = ^{ 2 , 3)5 however, ^^{ 2 , 3 } is nonempty only when X 2 < 0. Therefore, 
for a given departure point, depending on the coordinate X 2 , only one of these two 
subsets will be active with arrival points in different sectors Q{i, 2 , 3 } or Q{ 2 , 3 } (with 
different interpolation data). In a similar way, the remaining six subsets can be 
obtained. By assuming a linear control term, the identification of the control subsets 
is simple, as it will only require the resolution of linear inequalities where the departure 
point enters as a fixed data. 

Appendix B: Exact value function for Test 1. The exact value function 
for the infinite horizon optimal control problem with eikonal dynamics in Test 1 is 
derived. The HJB-equation has the form 

Av + max I-m^Vm- + =0, a; G M”, 7 > 0, 

ueu [ V 2 2 y I 


(4.4) 
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where U = {u G | \\u \\2 < 1}, for which we want to obtain an explicit solution. If 
“||V'r ’||2 < 1 , then = —^Vv provides a maximum for the expression in brackets, 
and the HJB-equation becomes 

Switching to polar coordinates (r, cp) this equation can be expressed as 
Xv{r, ip) + V>?) “ y = 0’ 

and assuming that the solution is radially symmetric 

1 

Xv{r) + - y = 0- 


The ansatz vi = Ar^, leads to XAr"^ + ^A'^r‘^ — ^ = 0 and hence A = 4 (-^A^ + ^ — A). 
The resulting expression for vi is the solution of ( |4.4| ), if {vi)r < 7 which results in 
r < 2((A^ + -)^ — A)“^ =: r. We note that vi{f) = 


We turn to the case that ^||V'i ;||2 > 1 . In this case the maximum in (4.4) is 


achieved on the boundary of U at 
metric solution in polar coordinates 


Vv 

■||v7|2 


Again we look for a radially sym- 


(4.5) 


Xv{r) + |n^(r)| “ y = |- 


Assuming that > 0 we make the Ansatz V 2 {r) = ar‘^ br cde Inserting 


into (4.5) and comparing coefficients we obtain 


/\ 7 f 7 —Xr 

Continuous concatenation with vi at r implies that 
(4.6) 


d = e 


Xr 


[(i-^ 

1 ^ 

\ - 1—2 

7 

1 ■ 



2A 

V. 


Note that this latter coupling yields that v{r) is a C^(M) function and therefore it is 
also a classical solution of eq. (4.4). Summarizing we have 


v{r) = 




7 2_7 

2a'’ A2'’ 


— -yr H-y + yy + de 


1 _ 

2 A 


1 

A3 


o-Ar 


for r < r 
for r > 


where d is given in (4.6). 
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Optimal trajectories and controls for initial position (0, l,7r/2) 



Fig. 4.5: Sparse control of a 3D car model. Top: different isosurfaces with different ^i-cost. 
Rows 2 and 3: effect of sparsity on the optimal control field. Bottom: optimal trajectories 
with and without sparsity. 

















